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Abstract. Perhaps the most extreme examples of "Active OB stars" are the subset of high-mass 
X-ray binaries - consisting of an OB star plus compact companion - that have recently been 
observed by Fermi and ground-based Cerenkov telescopes like HESS to be sources of very high 
energy (VHE; up to 30 TeV!) 7-rays. This paper focuses on the prominent 7-ray source, LS5039, 
which consists of a massive 06.5V star in a 3.9-day-period, mildly elliptical (e ~ 0.24) orbit 
with its companion, assumed here to be a black-hole or unmagnetized neutron star. Using 3-D 
SPH simulations of the Bondi-Hoyle accretion of the O-star wind onto the companion, we find 
that the orbital phase variation of the accretion follows very closely the simple Bondi-Hoyle- 
Lyttleton (BHL) rate for the local radius and wind speed. Moreover, a simple model, wherein 
intrinsic emission of 7-rays is assumed to track this accretion rate, reproduces quite well Fermi 
observations of the phase variation of 7-rays in the energy range 0.1-10 GeV. However for 
the VHE (0.1-30 TeV) radiation observed by the HESS Cerenkov telescope, it is important 
to account also for photon- photon interactions between the 7-rays and the stellar optical/UV 
radiation, which effectively attenuates much of the strong emission near periastron. When this is 
included, we find that this simple BHL accretion model also quite naturally fits the HESS light 
curve, thus making it a strong alternative to the pulsar-wind-shock models commonly invoked 
to explain such VHE 7-ray emission in massive-star binaries. 



1. Introduction 

Among the most active OB stars are those in high-mass X-ray binary (HMXB) systems, 
in which interaction of the massive-star wind with a close compact companion - either a 
neutron star or black hole - produces hard (> 10 keV) X-ray emission with characteristic, 
regular modulation over the orbital period. In recent years a small subset of such HMXB's 
have been found also to be gamma-ray sources, with energies up to ~10 GeV observed 
by orbiting gamma-ray observatories like Fermi, and very-high-energies (VHE) of over a 
TeV (10 12 eV) seen by ground-based Cerenkov telescopes like HESS, Veritas, and Magic. 
For the one case, known as B 1259-63, in which oberved radio pulses show the companion 
to be a pulsar, the 7— ray emission seems best explained by a Pulsar- Wind- Shock (PWS) 
model, wherein the interaction of the relativistic pulsar wind with the dense wind of 
the massive-star produces strong shocks that accelerate electrons to very high energies, 
with inverse- Compton scattering of the stellar light by these high-energy electrons then 
producing VHE 7-rays. For the other gamma-ray binaries, pulses have not been detected, 
and the nature of the companion, and the applicability of the PWS model, are less clear. 
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An altem&tiveMicroQuasar (MQ) model instead posits that accretion of circumstellar 
and/or wind material from the massive star onto the companion - assumed now to be 
either a black hole or a weakly magnetized neutron star - powers a jet of either relativistic 
protons, which interact with stellar wind protons to produce pions that quickly decay 
into 7-rays, or relativistic pair plasma, which inverse Compton scatters stellar photons 
to 7-rays. A poster paper by Okazaki et al. in these proceedings uses Smoothed Particle 
Hydrodyanmics (SPH) simulations to examine both the PWS and MQ models, applying 
them respectively for two systems, B 1259-63 and LSI +61 303, in which the massive star 
is a Be star. In such Be binary systems, the compact companion can interact with the Be 
star's low-density polar wind, and its dense equatorial decretion disk, and so modeling 
both types of interaction is key to to understanding their high-energy emission. 

In the paper here, we apply the same SPH code to a MQ model for the TeV-binary 
LS 5039, for which the (non-pulsing) compact companion is in a moderately eccentric 
(e ~ 0.24), 3.9-day orbit around a massive, non-Be primary star o f spectral type 06.5V . 
(See Table 1 for full parameters.) This builds on our previous study ( Okazaki et al.f 2008b) 



to account now for a fixed wind acceleration, instead of just assuming a constant wind 
speed. A key result is that the orbital variation of the 3-D SPH accretion rate of the stellar 
wind flow onto the compact companion follows very closely the analytic Bondi-Hoyle- 
Lyttleton (BHL) rate that depends on the wind speed and orbital separation at each 
phase. (See right panel of Figure [TJ) Assuming the mass accretion translates promptly 
into jet power and thus 7-ray emission, we then apply this result to derive predicted 
light curves at both the GeV energies oberved by Fermi and the TeV energies observed 
by HESS. For the former, we find that assuming GeV 7-ray emission track s closely the 



BHL accretion rate gives directly a quite good fit to the Fermi lightcurve ([Abdo et al 
2009. see left panel of Figure [3] below). But for the latter case one must also account 
for the attenuation of emitted TeV 7-rays by photon-photon interaction with the optical 
and UV radiation of the massive star. When this is in cluded, then the predi cted TeV 
lightcurve also closely matches the HESS observations ( Aharonian et al. 20061 see right 
panel of Figure [3] below). 

Table 1. Model parameters for LS 5039 

Primary Secondary 

Spectral type UB.5V Compact 

Mass (M Q ) 22. 9 a 3.7 a 

Radius 9.37?© a (= 0.31a) 2.5 x 10" 3 a 

Effective temperature T eB (K) 39,000 a 

Wind terminal speed Voo (kms -1 ) 2440 

Mass loss rate M* (M Q yr _1 ) 5 x 10" 7 a 

Orbital period P orb (d) 3.9060 a 

Orbital eccentricity e 0.24 b 

Semi-major axis a (cm) 2.17 x 10 12 

a Casares et al. ('2005) b Szalai et al. (2010) 



2. Bondi-Hoyle-Lyttleton (BHL) Accretion 

To provide a physical context for modeling accretion in LS 5039, it is helpful first 
to review briefly the basic scalings for accretion of an incoming flow onto a grav itating 
body , as analyzed in pioneering studies by Bondi, Hoyle and Lyttleton (BHL) (see Edgarl 



2004 for references and a modern review). The left panel of Figure [T] illustrates the basic 
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model. An initially laminar flow with density p and relative speed V re i is focussed by the 
gravity of a mass M onto the downstream side of the flow symmetry axis, whereupon 
the velocity component normal to the axis is cancelled by collisions among the streams 
from different azimuths. Material with initial impact parameter equal to a critical value 
b = 2GM/V^ el - now typically dubbed the BHL radius, and highlighted by the (red) 
heavy curve in Figure Q]- arrives on the axis a distance b from the mass, with a parallel 
flow energy V 2 e[ /2 just equal to the gravitational binding energy GM/b. Since all material 
with an impact < b (i.e. within the shaded area) should thus have negative total axial 
energy, it should be eventually accreted onto the central object. This leads to a simple 
BHL formula for the expected mass accretion rate, 



MbHL = pV re lTTb 2 = 



4tt P G 2 M 2 
K 3 el 
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Figure 1. Left: Illustration of flow streams for simple planar form of Bondi-Hoyle accretion. The 
gravitational attraction of mass M focuses an inflow with relative speed V re i onto the horizontal 
symmetry axis. This causes material impacting within a BHL radius b = 2GM /L r 2 ei of the axis 
to be accreted onto M, as outlined here by the bold (red) curve that bounds the lightly shaded 
region representing accreting material. Right: Orbital phase variation of the mass-accretion rate 
from 3-D SPH simulations (blue), compared to the analytic BHL formula from eqn. (|2.2|l (red) 
assuming the parameters given in Table [1] 



In a binary system like LS 5039 there arise additional effects from orbital motion and 
the associated coriolis terms. But if we ignore these and other complexities, we can use 
eqn. (|2.1[) to estimate how the accretion rate should change as a result of the changes in 
binary separation r over the system's elliptical orbit, 
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where M w = Airr 2 p(r)V w (r) is the wind mass loss rate, and V re i(r) is the local magnitude 
of the relative wind and orbital velocity vectors. We assume here that the stellar wind 
follows a standard 'beta-type' velocity law, V w (r) = V^o(l — i?»/r)' 3 , where R* is the 
O-star radius, is the wind terminal speed, and we adopt here a standard velocity 
power index (3 = 1. 

For the LS 5039 system parameters given in Table [I] the smooth (red) curve in the 
right panel of Figure [T] plots the orbital phase variation of this analytically predicted 
BHL-wind-accretion rate. A principal result of this paper is that this is in remarkably 
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close agreement with the jagged (blue) curve, representing the corresponding accretion 
rate from full, 3-D SPH simulations, the details of which we discuss next. 




-16 -15 -U -13 

log p 

Figure 2. Snapshots near periastron in SPH simulation of LS 5039, with the grayscale showing 
log of the density (in g/cm 3 ), and the arrows denoting flow direction and speed. The left panel 
gives an overall view of the O-star wind stream in the direction of the orbiting companion, which 
focuses the impingent wind flow into a narrow, dense, downstream wake. The right panel zooms 
in on this wake at a much smaller scale, within a BHL accretion radius b (marked by the dashed 
white circle arc) from the companion (represented by the black sphere of assumed accretion 
radius Race)- The white arrows show how dense material close to the companion and along the 
gravitational focus axis forms an accretion stream onto the companion, much as in the simple 
laminar flow picture illustrated in the left panel of Figure [T] 



3. SPH Simulations 

The SP H code used here is based on a version originally developed by Benz et all 



( 1990) and iBate et all (|l995h . with recent extensions to model interacting binaries by 



Okazaki et al. ( 2008al ). It uses a variable smoothing length and integrates the SPH equa- 



tions with the standard cubic-spline kernel using individual time steps for each particle. 
The artificial viscosity parameters are set to standard values of cusph = 1 and /?sph = 2. 
In the implementation here, the O-star wind is modeled by an ensemble of isothermal gas 
particles of negligible mass, while the compact object is represented by a sink particle 
of mass M and radius R aC c — 5 x 10 9 cm, i.e. much larger than the ~ 10 6 cm radius 
of a compact object, but still a factor 10 smaller than the minimum (periastron) value 
of the BHL accretion radius b; if SPH particles fall within this accretion sphere, they 
are removed from the simulation. The O star's mass exerts a gravitational pull on the 
binary companion, and the net force of radiative driving vs. gravity of the O-star leads 
to the outward wind acceleration characterized by the assumed beta=l velocity law. In 
addition, the individual SPH particles feel the gravitational attraction from the compact 
companion. Finally, to optimize the resolution and computational efficiency of our sim- 
ulations, the wind particles are ejected onl y in a narrow range of azimuthal and vertical 



angles toward the companion; figure 1 of lOkazaki et al.l (|2008bl ) shows that this gives 
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quite similar accretion flow structure to what is obtained in a full, spherically symmetric 
model, with however many fewer required particles. 

For the system parameters listed in Table Q] Figure [2] illustrates the nature of the 
accretion in these SPH simulations, using a time snapshot near periastron, with phase 
4> = 0.06. The left panel shows the overall, orbitally deflected wind stream that flows 
radially away from the O-star toward the companion, along with the narrow, dense, 
gravitationally focussed wind-stream in its wake. The right panel zooms in on this wake 
on a scale within a BHL accretion radius of the companion, as denoted by the dashed 
white circle arc. The white arrows show that, within the portions of this dense focal 
stream wake nearest the companion, the flow becomes directed toward the black sphere 
with assumed accretion radius Race- 

Averaged over some detailed variations, the overall process of accretion in this fully 3-D 
SPH simulation is thus remarkably similar to the simple laminar flow picture illustrated 
in the left panel of Figure [TJ The right panel of Figure Q] shows moreover that the SPH 
accretion rate - averaged over a phase interval 0.01 to smooth over rapid fluctuations - 
also agrees remarkably well with the simple BHL rate given by eqn. 
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Figure 3. LS 5039 light curves at energies 0.1-10 G eV observed by Fermi (|Abdo et al .1120091 . 
left), and at energies above 1 TeV observed by HESS ( Aharonian ct al. 200^, right), plotted vs. 
orbital phase, and compared against the model results based on BHL accretion, shown as heavy 
solid (blue) curves. (The lighter (red) curve represents the best-fit sinusoid for the HESS data.) 



4. Accretion-powered 7-ray emission 

In a microquasar model, the accretion onto the compact companion powers a high- 
energy jet that produces 7-rays. Since the jet acceleration occurs on a scale of a few tens 
of compact companion radii, it seems reasonable to postulate that the 7-ray emission 
could promptly track the accretion rate. The left panel of Figure [3] compares the Fermi 
lightcurve for LS 5039 at energies 0.1-10 GeV with a simple emission model that scales 
directly with the BHL accretion rate. The model fits the data quite well, reproducing 
both the roughly factor five variation range of observed 7-rays, as well as the observed 
flux peak at a phase very near periastroryj. 

By contrast, the right panel of Figure shows that the TeV 7-rays observed from 
LS 5039 by HESS peak well before periastron, toward the phase associated with inferior 
conjunction, with the flux minimum occuring just after periastron, near the phase for 
superior conjunction. These phase shifts in the peak and minimum can be explained by 

t If instead of a n eccentricity e = 0.24, we use the larger value e = 0.35 suggested by 
Casarcs ct al. (2005]), the amplitude of orbital variation in mass accretion rate then becomes 
stronger than this factor five 7-ray variation observed by Fermi. 
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accounting for the attenuation of the 7-ray rays through photon-photon interaction with 
the radiation from the O-star. For TeV 7-rays, the geometric mean of the energies of the 7- 
rays and stellar UV/optical photons exceeds twice the rest-mass-energy of electrons, thus 
allowing photon-photon production of electron/positron pairs. The thick (blue) curve in 
Figure [3] shows that a simple model combining a BHL-rate emission with such photon- 
photon attenuation does indeed match the HESS data very well. In contrast, for GeV 
energy 7-rays the geometric mean with stellar UV/optical light generally falls below the 
electron mass-energy, and so the Fermi lightcurve is unaffected by such attentuation. 

Since the cross section for photon-photon interaction is highest near the threshhold of 
about 0.1 TeV, a pure attenuation model predicts that the HESS spectrum should be 
hardest during the superior conjunction phase of minimum flux and maximum attenu- 
ation. However, plots of the photon index for TeV observations by HESS show just the 
opposite, with the minimum flux corresponding to the softest spectrum. We are currently 
investigating whether this can be explained by accounting for the progressive softening 
of 7-rays during a cascade of absorption and reemission associated with pair production 
and anniliation. 
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Discussion 

DOUG Gies: In Bondi-Hoyle accretion, the gas falls towards the black hole and presum- 
ably forms a small accretion disk that powers the relativistic jets. In your model, how 
long does the gas reside in the disk before being launched into the jets? 

Stan Owocki: Our modeling thus far effectively assumes this residence time is short, 
but this requires further investigation through a more detailed model of the accretion 
powering of the jet. 

Guillaume DuBUS: I know of no accretion model showing that gamma ray emission 
from accretion can promptly follow the orbital variation of the BHL accretion rate. Also, 
how can the MQ model explain the low fraction of gamma-ray sources among X-ray 
binaries? The PWS model explains this through the fact they exist for only a short time 
in a state that can give the required relativistic wind. I thus believe such PWS models 
are much favored over an MQ model. 
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Stan Owocki: Our co-author on this paper, Gustavo Romero, who could not attend 
this meeting, has voiced to us similarly strong arguments in favor of the MQ model over 
the PWS model for LS5039. I expect there will be a vigorous debate on these issues at 
the upcoming gamma-ray meeting in Heidelberg this December. 



